/*这里我们进行多种计算，我们知道在cuSPARSE库中，主要有三个级别的乘法计算
  https://www.jianshu.com/p/0c7e8848fd11
  一个级别我们举的例子是将稀疏向量x中的元素散步到向量y中，并且覆盖y中已有的向量

*/
#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<cusparse.h>

int main(){
  //我们首先在CPU上创建一个稀疏矩阵和一个线性矩阵
  /*
  x = [100 , 200 , 0 , 400]，我们就装它是一个稀疏矩阵
  y = [1 , 2 , 3 , 4 , 5 , 6 , 7 , 8]，这是一个稠密矩阵
  我们要将x中的元素散发到&y[4]中，使用接口cusparse<t>sctr()
  这样子会覆盖y中的元素
  最终结果y = [1 , 2 , 3 , 4 , 100 , 200 , 7 , 400]
  */

  //我们创建一个数组
  float* yHost = 0;
  float* xHostVal = 0;
  int* xHostCol = 0;
  //y向量的长度
  int ySize = 8;
  //x中非零元素的数量
  int xnzSize = 3;

  //我们创建y数组
  cudaMallocHost((void **)&yHost, ySize * sizeof(float));
  
  //用来遍历的变量
  int i = 0;
  for(i = 0; i < ySize; i++){
    yHost[i] = i + 1;
  }

  //我们创建x数组，我们使用压缩存储，一个存储x中的值，一个是x的非零元素所在列
  cudaMallocHost((void **)&xHostVal, xnzSize * sizeof(float));
  cudaMallocHost((void **)&xHostCol, xnzSize * sizeof(int));

  xHostVal[0] = 100;
  xHostCol[0] = 0;
  xHostVal[1] = 200;
  xHostCol[1] = 1;
  xHostVal[2] = 400;
  xHostCol[2] = 3;

  //然后我们将这些数组传到设备端
  //创建设备端数组
  float* yDevice = 0;
  float* xDeviceVal = 0;
  int* xDeviceCol = 0;
  cudaMalloc((void **)&yDevice, ySize * sizeof(float));
  cudaMalloc((void **)&xDeviceVal, xnzSize * sizeof(float));
  cudaMalloc((void **)&xDeviceCol, xnzSize * sizeof(int));

  //将数据从主机端传到设备端
  cudaMemcpy(yDevice, yHost, ySize * sizeof(float), cudaMemcpyHostToDevice);
  cudaMemcpy(xDeviceVal, xHostVal, xnzSize * sizeof(float), cudaMemcpyHostToDevice);
  cudaMemcpy(xDeviceCol, xHostCol, xnzSize * sizeof(int), cudaMemcpyHostToDevice);

  //CUDA官方库的使用方法和面向对象非常像，我们首先要创造一个cuSPARSE的句柄
  cusparseHandle_t handle = 0;
  cusparseCreate(&handle);

  //然后使用cuda的库函数来执行稀疏矩阵x非零元素在稠密矩阵y中的散发操作，cusparseSsctr
  //CUSPARSE_INDEX_BASE_ZERO代表矩阵的索引从0开始记录
  cusparseSsctr(handle, xnzSize, xDeviceVal, xDeviceCol, &yDevice[ySize/2], CUSPARSE_INDEX_BASE_ZERO);
  cudaThreadSynchronize();
  //将数据传回
  cudaMemcpy(yHost, yDevice, ySize * sizeof(float), cudaMemcpyDeviceToHost);

  for(i = 0; i < ySize; i++){
    printf("%f,", yHost[i]);
  }

  printf("\n");

  //回收主机端空间
  cudaFreeHost(yHost);
  cudaFreeHost(xHostVal);
  cudaFreeHost(xHostCol);

  //回收设备端空间
  cudaFree(yDevice);
  cudaFree(xDeviceVal);
  cudaFree(xDeviceCol);
}